home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / SPLINES / EX / SPLINEIN.C < prev   
Encoding:
C/C++ Source or Header  |  1996-04-18  |  4.5 KB  |  125 lines
  1. #include <BSpline.h>
  2. #include <ElmItgRules.h>
  3. #include <math.h>
  4.  
  5. // This demo program is meant to demonstrate some features of
  6. // the spline package and class 'ElmItgRules' in DIFFPACK.
  7.  
  8. // Two functions (sin(x) & cos(x) for ex.) are interpolated by splines
  9. // 'spline1' & 'spline2', then \int_a^b (spline1*spline2) is calculated.
  10. // Numerical integration is used in the calculation.
  11. // The optimization due to the unform knot vector is not implemented.
  12.  
  13. // ======================== optimal version ==========================
  14.  
  15. // specific constants for the particular case of the f_1 and f_2 functions:
  16. const real a = 0.0;          // integration region
  17. const real b = M_PI/2.0;
  18. int k = 4;                   // cubic splines as default
  19. int nitg_order = 4;          // order of numerical integration
  20.  
  21. real f_1 (real x) { return sin(x); }  // two analytical functions
  22. real f_2 (real x) { return cos(x); }
  23.  
  24. int main (int nargs, const char** args)
  25. {
  26.   if (nargs<2)
  27.     errorFP("main:","Usage: %s num_knots (nitg_order k)",args[0]);
  28.  
  29.   if (nargs>2) nitg_order = atoi(args[2]);
  30.   if (nargs>3) k = atoi(args[3]);
  31.  
  32.   KnotVec knots (atoi(args[1])); // number of unique knots from input
  33.   knots.fill (a,b);              // uniform knot vector (not necessary)
  34.   knots.regulate (k);            // make a k-regular knot vector
  35.  
  36.   SplineSpace space (knots,k);   // the corresponding spline space
  37.  
  38.   int n = space.dimension();     // total number of B-spline basis functions
  39.   int derivative_order_i=0;      // derivative order for B_i(x)
  40.   int derivative_order_j=0;      // derivative order for B_j(x)
  41.   Mat(real) BS(1,k);             // assistant variables
  42.   int miu1,miu2,l,r,i,j,m,s;
  43.   real x;                        // a point in [a,b]
  44.   real scale;                    // scaling for x -> [-1,1] num.itg. interval
  45.  
  46.   BSpline spline1, spline2;      // two spline functions
  47.   spline1.redim (space);
  48.   spline2.redim (space);
  49.  
  50.   // let spline1 interpolate the f_1 function and spline2 interpolate f_2:
  51.   Vec(real) xs(n);               // x coordinates in [a,b]
  52.   Vec(real) fs(n);               // f_1 and f_2 values over xs
  53.   Vec(real) coef1(n), coef2(n);  // spline coeff. for spline1 and spline2
  54.   xs.fill (a,b);
  55.   fs = xs; fs.apply (f_1);       // fs contains f_1 values
  56.   spline1.interpolation (xs,fs); // spline interpolation of function f_1
  57.   spline1.getCoef (coef1);       // coefficients for spline1
  58.   fs = xs; fs.apply (f_2);       // fs contains f_2 values
  59.   spline2.interpolation (xs,fs); // spline interpolation of function f_2
  60.   spline2.getCoef (coef2);       // coefficients for spline2
  61.  
  62.   ElmItgRules rule (GAUSS_POINTS); // used for numerical integration
  63.   rule.refill (nitg_order, 1 /*1D*/, BOX /*domain: box (or interval)*/);
  64.   Ptv(real) pt(1);               // used to hold a numerical integration point
  65.  
  66.   // B-spline values at all the numerical integration points:
  67.   ArrayGenSimple(real) bs1(n,nitg_order,k), bs2(n,nitg_order,k);
  68.   bs1.fill(0.0); bs2.fill(0.0);
  69.  
  70.   // first calculate all B-spline values, bs1 and bs2:
  71.   for (m=1; m<=n; m++)
  72.     if (knots(m)<knots(m+1))   // nonempty interval
  73.       {
  74.         miu1 = miu2 = m;
  75.         scale = (knots(m+1)-knots(m))/2.0;
  76.  
  77.         for (s=1; s<=nitg_order; s++)
  78.           {
  79.             rule.point (pt,s);
  80.             x = knots(m) + (pt(1)+1.0)*scale;
  81.  
  82.             space.BSplines (x,derivative_order_i,BS,miu1,dpTRUE);
  83.  
  84.             for (i=1; i<=k; i++) bs1(m,s,i) = BS(1,i);
  85.  
  86.             if (derivative_order_j != derivative_order_i)
  87.               space.BSplines (x,derivative_order_j,BS,miu2,dpTRUE);
  88.  
  89.             for (i=1; i<=k; i++) bs2(m,s,i) = BS(1,i);
  90.           }
  91.       }
  92.   real I = 0.0;              // for containing integration result
  93.  
  94.   for (i=1; i<=n; i++)       // integrate numerically
  95.     for (j=1; j<=((derivative_order_i == derivative_order_j) ? i:n); j++)
  96.       if ( abs(i-j) < k )    // B_i(x)*B_j(x) is different from 0
  97.         {
  98.           l = max(i,j);      // actual integration region
  99.           r = min(i,j)+k;
  100.  
  101.           for (m=l; m<r; m++)
  102.             if (knots(m+1)>knots(m))
  103.           {
  104.         scale = (knots(m+1)-knots(m))/2.0;
  105.         miu1 = miu2 = m;
  106.  
  107.         for (s=1; s<=nitg_order; s++)
  108.           {
  109.             if (derivative_order_i==derivative_order_j && j < i)
  110.               I += (coef1(i)*coef2(j)+coef1(j)*coef2(i))
  111.             *bs1(m,s,i-miu1+k)*bs2(m,s,j-miu2+k)
  112.               *rule.weight(s)*scale;
  113.             else
  114.               I += coef1(i)*coef2(j)
  115.             *bs1(m,s,i-miu1+k)*bs2(m,s,j-miu2+k)
  116.               *rule.weight(s)*scale;
  117.           }
  118.           }
  119.         }
  120.  
  121.   s_o<<oform("Integration result with spline k=%d and nitg_order=%d is %g\n",
  122.              k,nitg_order,I);
  123.   return 0;
  124. }
  125.